arXiv: 1505.00792 v 3 [astro-ph.HE] 21 Dec 2015 


Binary Black Hole Mergers from Globular Clusters: Implications for Advanced LIGO 


Carl L. Rodriguez, 1 Meagan Morscher, 1 Bharath Pattabiraman, 1 ’ 2 
Sourav Chatterjee, 1 Carl-Johan Haster, 1,3 and Frederic A. Rasio 1 

1 Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) and Dept, of Physics and Astronomy, 
Northwestern University, 2145 Sheridan Rd, Evanston, IL 60208, USA 
2 Dept, of Electrical Engineering and Computer Science, Northwestern University, Evanston, IL, USA 
3 School of Physics and Astronomy, University of Birmingham, Birmingham, B15 2TT, United Kingdom 

(Dated: December 22, 2015) 

The predicted rate of binary black hole mergers from galactic fields can vary over several orders 
of magnitude and is extremely sensitive to the assumptions of stellar evolution. But in dense stellar 
environments such as globular clusters, binary black holes form by well-understood gravitational 
interactions. In this letter, we study the formation of black hole binaries in an extensive collection 
of realistic globular cluster models. By comparing these models to observed Milky Way and extra- 
galactic globular clusters, we find that the mergers of dynamically-formed binaries could be detected 
at a rate of ~ 100 per year, potentially dominating the binary black hole merger rate. We also find 
that a majority of cluster-formed binaries are more massive than their field-formed counterparts, 
suggesting that Advanced LIGO could identify certain binaries as originating from dense stellar 
environments. 


INTRODUCTION 

By the end of this decade, the Advanced LIGO and 
Virgo detectors are expected to observe gravitational 
waves (GWs), ushering in a new post-electromagnetic 
era of astrophysics mm- The most anticipated sources 
of observable GWs will be the signals generated by merg¬ 
ers of binaries with compact object components, such as 
binary neutron stars (NSs) or binary black holes (BHs). 
While coalescence rates of NS-NS or BH-NS systems can 
be constrained from observations, it is not currently pos¬ 
sible to produce observationally-motivated rate predic¬ 
tions for BH-BH mergers [3J. Typical detection rates of 
binary BH (BBH) mergers in galaxies can span several 
orders of magnitude from 0.4 yr _1 to 1000 yr _1 with 
a fiducial value of ~ 20 yr _1 0; however, these esti¬ 
mates typically ignore the large numbers of BBHs that 
are formed through dynamical interactions in dense star 
clusters 00. 

The dynamical formation of BBHs is a probabilistic 
process, requiring a very high stellar density. These con¬ 
ditions are believed to exist within the cores of globular 
clusters (GCs), very old systems of ~ 10 5 —10 6 stars with 
radii of a few parsecs. Approximately 10 Myr after the 
formation of a GC, the most massive stars explode as su¬ 
pernovae, forming a population of single and binary BHs 
with individual masses from ~ 5 Mq to ~ 25 Mq [7]. The 
BHs, being more massive than the average star in the 
cluster, sink to the center of the GC via dynamical fric¬ 
tion, until the majority of the BHs reside in the cluster 
core 0. After this “mass segregation” is complete, the 
core becomes sufficiently dense that three-body encoun¬ 
ters can frequently occur 0, producing BBHs at high 
rates. In effect, GCs are dynamical factories for BBHs: 
producing large numbers of binaries within their cores 
and ejecting them via energetic dynamical encounters. 


In this letter, we use an extensive and diverse collection 
of GC models to study the population of BBHs that Ad¬ 
vanced LIGO can detect from GCs. We explore how the 
observed parameters of a present-day GC correlate with 
the distribution of BBH inspirals it has produced over its 
lifetime. We then compare our models to the observed 
population of Milky Way GCs (MWGCs) and use recent 
measurements of the GC luminosity function to deter¬ 
mine a mean number of BBH inspirals per GC. Finally, 
we combine these estimates with an updated estimate of 
the spatial density of GCs in the local universe (Appendix 
I) into a double integral over comoving volume and in¬ 
spiral masses to compute the expected Advanced LIGO 
detection rate. We assume cosmological parameters of 
Om = 0.309, Ha = 0.691, and h = 0.677, consistent with 
the latest combined Planck results ljjj. 


COMPUTING THE RATE 

We use a collection of 48 GC models generated by 
our Cluster Monte Carlo (CMC) code, an orbit-averaged 
Henon-type Monte Carlo code for collisional stellar dy¬ 
namics El- The models span a range of initial star 
numbers (2 x 10 5 to 1.6 x 10 6 ), initial virial radii (0.5 
pc to 4 pc), and consider low stellar metallicities (Z = 
0.0005, 0.0001) and high stellar metallicities {Z = 0.005). 
In addition, the code implements dynamical binary for¬ 
mation via three-body encounters, strong three and four- 
body binary interactions, and realistic single and binary 
stellar evolution. See Appendix II for a complete descrip¬ 
tion of our code and the models used. 

Previous studies have explored the contribution of 
BBHs from GCs to the Advaned LIGO detection rate 
HMZj; however, the majority of these studies have re¬ 
lied on either approximate analytic arguments or sirnpli- 
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fied iV-body models with N < 10 s particles and have 
assumed a single black hole mass of IOMq. The one ex¬ 
ception is [15], which used a Monte Carlo approach to 
model GCs of a realistic size (N = 5 x 10 5 ). However, 
their study only considered GCs of a single mass, and 
did not extrapolate that result to the observed distribu¬ 
tion of GCs in the local universe. Ours is the first study 
to compare models with all the relevant physics over a 
range of masses to observed GCs. This comparison en¬ 
hances our BBH merger rate by more than an order of 
magnitude over previous results. 

We express the rate of detectable mergers per year, Rd, 
as the following double integral over binary chirp mass 
(M c = (mim 2 ) 3/5 /(mi + TO 2 ) 1 / 5 ) and redshift: 

Rd = JJ U(M c ,z)f d (M c ,z) < ^ ( ^dM c dz . (1) 

This equation is similar to that found in mm- The 
components of Eqn.[l]are as follows: 

• 1Z(A4 C , z ) is the rate of merging BBHs from GCs 
with chirp mass A4 C at redshift z. 

• z ) is the fraction of sources with chirp mass 
M c at redshift z that are detectable by a single 
Advanced LIGO detector. 

• dV c /dz is the comoving volume at a given redshift 

m- 

• dtg/dto = 1/(1 + z) is the time dilation between 
a clock measuring the merger rate at the source 
versus a clock on Earth. 

This letter focuses on estimating the rate, TZ(Ai c ,z), 
using our collection of GC models. We assume the rate 
can be expressed as the product of the mean number 
of inspirals per GC, the distribution of those sources in 
M c — z space, and the density of GCs in the local uni¬ 
verse, i.e. lZ(M c ,z) = ( N ) x P(M c ,z) x pac- The spa¬ 
tial density of GCs in the local universe is taken to be 
Pgc = 0.77 Mpc -3 , based on recent measurements of 
extragalactic GC systems [21] and modern near-infrared 
Schechter functions [22] • Note that this estimate, com¬ 
puted in Appendix I, is substantially lower than the pre¬ 
vious estimate of pac = 8.4 ft 3 Mpc -3 from [T2] that has 
been used in previous studies. We now estimate the val¬ 
ues of ( N ) and P(A4 c ,z). 


Inspirals vs GC Mass 



FIG. 1. The number of BBH inspirals per model at 12 Gyr for 
each of our 48 GC models as a function of final cluster mass. 
We show the weighted linear regression (with la uncertainties 
on the slope) for the low and high metallicity models. 

we compare the total masses and concentrations of our 
models to GCs observed in the Milky Way. The con¬ 
centrations are measured by considering the ratio of a 
cluster’s core radius to its half-light radius, R c /Rh ■ This 
mass-concentration space is similar to the “fundamental 
plane” of GCs described in [23], with R c /Rh in place of 
the King concentration [23] , 

Two trends emerge in our models. First, the total 
number of BBH inspirals over 12 Gyr is nearly linearly 
proportional to the final cluster mass. Second, the num¬ 
ber of inspirals is higher for more compact clusters (those 
with smaller R c /Rh). Since the model coverage of the 
R c /Rh space is poorer than the coverage of the mass, and 
since there are no detailed observations of extragalactic 
GC concentrations, we elect to focus on the linear re¬ 
lationship between a GC’s mass and the number of in¬ 
spirals it has produced. We perform a weighted linear 
regression for both low-metallicity and high-metallicity 
GCs (Fig. [l]). The weights are created by generating a 
kernel density estimate (KDE) [25] of the MWGCs in 
the fundamental plane, then measuring the probability 
of each model as reported by the KDE. In other words, 
GC models that are more likely to represent draws from 
the distribution of MWGCs are more heavily weighted. 

We compute the mean number of inspirals per GC by 
multiplying the linear relationships from Fig. ]T] by the 
mass distribution of GCs. Recent work [26] has suggested 
that the distribution of GC luminosities is universal and 
well-described by a log-normal distribution: 


MEAN NUMBER OF MERGERS PER CLUSTER 

To determine the mean number of BBH inspirals pro¬ 
duced by a GC, we use the collection of models to ex¬ 
plore how the present-day observable parameters of GCs 
relate to the number of BBHs it has produced over its 
lifetime. To quantify the realism of a particlar model, 


dN 

d\ogL 


N 0 exp - 


(log T — log£ 0 )" 
2cr£ 


( 2 ) 


with logLo = 5.24 and = 0.52. Assuming a mass- 
to-light ratio of 2 in solar units 123122, we convert this 
luminosity function to a mass function. We then compute 
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the average number of inspirals per GC by integrating the 
linear relationship over the normalized GC mass function. 
The results for both metallicities and different high-mass 
cutoffs are shown in Table U) 


MWGCs will have large numerators and large denomina¬ 
tors, yielding smaller weights; however, as we will select 
some number of inspirals from each of these neighboring 
models, the cumulative effect is the same. 



Mass Cutoff 

Metallicity 

4 x 1O 6 AL 0 

2 x 1O 7 M 0 

2 x 1O 8 M 0 

Low 

430 

967 

1512 

High 

830 

1954 

3103 


TABLE I. The mean number of inspirals per GC over 12 Gyr. 
The result depends on our choice of maximum GC mass. We 
consider cutoffs of 4 x 1O 6 M0 (the approximate mass of the 
most massive MWGC, w Cen), 2 x 10 7 Mq (the approximate 
cutoff used in [26]), an d 2 x 1O 8 M 0 (the mass of the ultra¬ 
compact dwarf M60-UCD1 [28)1. 


DISTRIBUTION OF INSPIRALS 

The numbers quoted in Table [T] provide us with the 
mean number of BBH inspirals from a GC over 12 Gyr. 
We could use this average rate to compute a detection 
rate for Advanced LIGO. However, it is qualitatively ob¬ 
vious that the mass distribution of BBH sources is not 
constant in time (Fig. [2]). 

Therefore, we must use the distribution of BBH in¬ 
spiral events over time from GCs to compute the rate. 
We select inspirals randomly from each of our models, 
drawing more inspirals from models with higher weights 
according to the following scheme: 


W(M, R c /R h ) 


KDE M \ngc(M, Rc/R h) 

K DE Mo deis{M , RJ R h ) 


(3) 


where the weight, W{M,R C /Rh), of a model with mass 
M and compactness R c /Rh at 12 Gyr is defined by the 
ratio of the MWGC KDE at M,R C /Rh divided by the 
KDE of the models themselves, evaluated at M,R C /Rh- 
The reason for these weights is as follows: we wish to 
draw more samples from models that are more likely to 
represent MWGCs, but because our collection of models 
is drawn from a different distribution (the initial condi¬ 
tions from 0 ), we cannot simply draw inspirals at ran¬ 
dom from each model according to how well it repre¬ 
sents real GCs. To do so would bias our samples with 
the distribution that results from our initial conditions. 
By dividing the probability of a model representing a 
MWGC by the probability density of our collection of 
models, our scheme naturally corrects for this. Mod¬ 
els unlikely to represent MWGCs have small numerators 
and low weights. Models with no neighboring models 
that are likely to represent MWGCs have large numera¬ 
tors and small denominators, yielding high weights. Con¬ 
versely, models with neighbors that are likely to represent 
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Redshift 


FIG. 2. A sampe distribution of inspirals in redshift from 
the set of models. The redshift is computed by assuming 
that the difference between the present day and the inspiral 
time corresponds to the cosmological lookback time at a given 
redshift (e.g. 20). The number of inspirals drawn from each 
model is proportional to its weight, or how similar it is to 
the observed distribution of MWGCs. Inspirals of BBHs that 
were formed primordially are indicated with stars (merged in 
the cluster) and diamonds (ejected before merger). Inspirals 
of BBHs formed dynamically are shown as squares (in-cluster) 
and circles (ejected). Note that there are no binaries that are 
formed by binary stellar evolution with chirp masses greater 
than ~ 13M 0 (dashed line). This result is consistent across 
all models. The blue shaded regions illustrate the regions 
of parameter space where 50%, 10%, and 1% of sources are 
detectible by Advanced LIGO. 


We show a sample distribution of the chirp masses 
versus redshift in Fig. [2] We distinguish between two 
different BBH formation channels: primordial and dy¬ 
namical. We define primordial BBHs as those that are 
formed from the supernovae of two main sequence stars 
in a binary, and whose components were never bound 
to any other star before merger; conversely, we define 
dynamical binaries as those that are either formed from 
two isolated BHs via a three-body encounter, or formed 
from a higher-order dynamical encounter (a binary-single 
or binary-binary interaction forming a new binary pair). 
Primordial binaries can still have their orbital parame¬ 
ters modified by dynamics (via a strong encounter with 
another BH or BBH), as long as the encounter leaves 
the primordial BBH intact. One immediately appar¬ 
ent feature is the bi-modality between primordial and 
dynamical BBHs. Over all of our models, the highest 
chirp mass that is formed by pure binary stellar evolu¬ 
tion is Ai c ~ 13-Mq, as systems with larger progenitors 
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are distrupted by the supernova kick. This implies that 
any source from our models with a detected chirp mass 
greater than ~ 13M 0 could only have formed dynami¬ 
cally. 

To compare this result to BBHs formed in the field, we 
generated two additional CMC models, each containing 
5 x 10 6 binaries and different metallicites (Z = 0.005 and 
Z = 0.0005). These models were computed without two- 
body relaxation, binary formation, or strong encounters, 
and only considered the physics of binary stellar evolu¬ 
tion. In this dynamics-free environment, the maximum 
chirp mass of any BBH was M. c ~ 13 Mg. Although this 
result depends on the metallicity and the physics of stel¬ 
lar evolution, it does suggest that GC dynamics forms 
BBHs consistently more massive than those in the field 

DETECTION RATE 

We now compute the expected rate of signals de¬ 
tectable by Advanced LIGO. To compute the fraction 
of detectable sources, fd(M c ,z), we use gravitational 
waveforms that cover the inspiral, merger, and ringdown 
phases of a compact binary merger (known as IMRPhe- 
nomC waveforms [29]) and compute the signal-to-noise 
ratio (SNR) using the projected zero-detuning, high- 
power configuration of Advanced LIGO [3D]. We then 
marginalize over binary orientation and sky location to 
determine what fraction of sources at a given chirp mass 
and redshift yield an SNR > 8. This approach is identi¬ 
cal to that found in m- Note that we have assumed all 
binary components have equal masses in order to sim¬ 
plify the integral. This assumption is well-justified, as 
the dynamically-formed BBHs in our models have sim¬ 
ilar component masses. We assume all BHs to be non¬ 
spinning. 

Our distribution of inspirals, P(A4 c ,z), is generated 
by creating a KDE of the inspirals in Fig. [2] Since each 
“draw” of inspirals will produce a slightly different distri¬ 
bution, we compute P(A4 c ,z) 1000 times, and then take 
the mean of Eqn. [T]for those 1000 draws. We find that a 
single Advanced LIGO detector operating at design sen¬ 
sitivity will detect ~ 100 BBHs per year from GCs. Of 
those about 2/3 will originate from low-metallicity GCs 
and the rest from high-metallicity GCs, assuming 76% of 
GCs are low metallicity (consistent with the MWGC dis¬ 
tribution). Approximately 80% of these sources will have 
chirp masses greater than 13(1 + -z)Mq, meaning that the 
majority of BBHs detectable by Advanced LIGO from 
GCs could only havie formed dynamically. 

The majority of these BBH sources will be detected at 
low redshifts. For low-metallicity clusters, the distribu¬ 
tion of detectible sources in redshift peaks at z ~ 0.3, 
while for high-metallicity clusters the distribution peaks 
at z ~ 0.24. In both cases, 90% of detectable sources are 
located at z < 0.57. 


To obtain a rough estimate of the uncertainty on this 
prediction, we perform a simple error analysis that con¬ 
siders the optimistic and conservative rates that would 
be obtained by varying our assumptions and selecting 
the drier estimates of certain quantities. For the conser¬ 
vative estimate, we assume that the GC mass function 
truncates at the mass of ui Cen (4 x 10 6 Mq), and that 
the spatial density of GCs is pec = 0.32 Mpc -3 (the con¬ 
servative estimate from Appendix I). We also recompute 
the rate using the — la uncertainty from the regression 
in Fig. |T] and the lower standard deviation of our 1000 
draws of lZ(M c ,z). This yields a conservative estimate 
of ~ 20 BBH inspirals per year. Conversely, if we as¬ 
sume the most optimistic truncation mass for GC mass 
function (2 x 10 8 Mg), the most optimistic GC spatial 
density (pac = 2.3 Mpc -3 , the optimistic estimate from 
Appendix I, similar to the value used in previous stud¬ 
ies), and the +la uncertainties on the linear regression 
and 1Z(M C , z), we find an optimistic rate of ~ 700 BBH 
inspirals per year. This range is primarily influenced by 
the uncertainty in the GC spatial density and the trun¬ 
cation mass of the GC mass function. 


CONCLUSION 

In this letter, we compared new GC models computed 
with our CMC code to the observed distributions of 
Milky Way and extragalactic GCs to predict the expected 
rate of BBH inspirals from realistic GCs. We determined 
a linear relationship between the present-day mass of a 
GC and the number of BBH inspirals produced by that 
cluster. By combining this with the universal GC lumi¬ 
nosity function and a new estimate for the spatial den¬ 
sity of GCs, we were able to predict the mean density 
of BBH inspirals from GCs in the local universe. Then 
by weighting our models according to their similarity to 
the observed distribution of MWGCs, we created a dis¬ 
tribution of inspiral sources in chirp mass and redshift. 
Finally, by combining this with the anticipated sensitiv¬ 
ity of Advanced LIGO, we estimated a detection rate of 
~ 100 BBH inspiral events per year from GC sources. 
With highly conservative assumptions, this rate drops to 
~ 20 events per year, while highly optimistic assumptions 
pushes the rate as high as ~ 700 events per year. 

We also found that no BBHs with chirp masses above 
~ 13 Mq were formed from a primordial binary. In other 
words, every inspiral with A4 C > 13Mg in our models 
was formed by dynamical processes alone. This could, in 
theory, provide an easy way to discriminate between bi¬ 
naries that were formed dynamically versus those formed 
by binary stellar evolution; however, this result is highly 
dependent on the physics of supernova kicks and the frac¬ 
tion of ejected supernova material which falls back onto 
the newly formed BH, both of which remain poorly con¬ 
strained. In addition, recent work has suggested that the 
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mass distribution of chirp masses for BBHs produced by 
stellar evolution can reach as high as Ai c ~ 30M©, de¬ 
pending on the physics of the common envelope |31j . As 
such, this result should be treated as a proof-of-principle, 
and not a concrete physical claim. Investigations to bet¬ 
ter understand the relationship between this formation 
cutoff, the distribution of supernovae kicks, and the fall¬ 
back fraction, are currently underway. 

Finally, the number of BHs formed is entirely depen¬ 
dent on the choice of the initial mass function. Although 
our choice of IMF is typical for studies of this type, a 
variation of lcr in the slope of the high-mass end of the 
IMF can produce significant differences in the number 
of BBHs. Investigations to quantify this effect are also 
underway. 
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Number of Globular Clusters per Mpc 3 

In order to estimate the rate of inspirals from an aver¬ 
age GC per Mpc 3 , we must compute the average spatial 
density of GCs in the local universe. This can be ac¬ 
complished by considering the mean number of GCs per 
galaxy at a given luminosity, multiplied by the spatial 
density of galaxies at that luminosity, then summing over 
all luminosities, as illustrated in the following equation: 


PGC 


f # of GCs \ /# of Galaxies \ r 

\Galaxy/Mag ) X \ Mpc 3 x Mag / &S ' 


(4) 

The number of GCs per galaxy per luminosity can be 
determined by use of the Harris Globular Cluster System 
catalog [21]. The catalog provides a list of 422 galax¬ 
ies, their morphological type, visual and K-band magni¬ 
tudes (where available), and the estimated total number 
of GCs. In Figure [3] we plot the 346 galaxies for which 
K-band photometry is available in the catalog against the 
estimated number of GCs. For each collection of galaxy 
morphologies, we perform a Gaussian Process regression 
with the George package, described in [31]. The Gaus¬ 
sian processes are generated using a squared-exponential 
kernel combined with a white noise kernel, and then fit to 
the log of the number of GCs per galaxy. The kernel hy¬ 
perparameters are selected by maximizing the marginal¬ 
ized log-likelilrood of the Gaussian Process. See [33] for 
a detailed description of regression with Gaussian Pro¬ 
cesses. The mean and standard deviation of the resulting 


Gaussian Processes are also shown in Fig. [3] Note that 
the catalog does not include low-luminosity dwarf and 
irregular galaxies for which Ngc = 0. This suggests that 
our fitted function can be systematically overestimating 
the number of GCs from dwarf early-type galaxies with 
M v > -18 [M|. 


nr 




icr 


icr 


10 


io u 



M k 



M k M k 


FIG. 3. Number of GCs per galaxy per K-band luminosity 
for different galaxy morphologies in the Harris GC Systems 
Catalog m- The top plot shows the fit to all galaxies in the 
catalog, and is the regression used for the spatial estimate of 
GCs in the text. The red line shows the mean of the Gaus¬ 
sian process regression for each sample, while the gray shows 
the Ict confidence interval about each estimate. Note that 
this catalog does not include observed dwarf galaxies with no 
GCs, suggesting that the estimated mean at low magnitudes 
is systematically biased to higher values. 


The mean in Figure [3] gives us the number of GCs 
per galaxy per K-band luminosity. To compute a spatial 
density, we must integrate the mean over the density of 
galaxies per K-band luminosity per Mpc -3 . For this, we 
can use the well known Schechter function [35] , which de¬ 
scribes the spatial density of galaxies per absolute magni¬ 
tude interval dM. We use the recently-computed K-band 
Schechter functions fits from the Galaxy and Mass As- 
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sembly (GAMA) Survey [22], which contains individual 
fits to individual Schechter functions for each galaxy mor¬ 
phology, and a single fit to a double Schechter function 
for the combined sample of all galaxies. We use both 
to determine our overall density of GCs, as well as the 
density contributed by galaxies of each type. 

Finally, we multiply each Schechter function from [22] 
by the estimated number of GCs per galaxy from Fig¬ 
ure [3| and integrate over all K-band magnitudes (M k < 
— 15). This results in a GC density estimate of p G c = 
0.77 Mpc -3 , which we employ in our rate calculation. For 
completeness, we also consider different cutoffs for our 
magnitude integral, and report the contribution to the 
spatial densities from each galaxy morphology, in Table 

m 

Since the Schechter Function diverges at low lumi¬ 
nosities, and since our fit systematically overestimates 
the number of GCs for low-luminosity galaxies, we must 
pick a reasonable limit at which to truncate our integral. 
We use a lower limit of M k = —15, although for com¬ 
pleteness we also consider lower (M k = — 17) and higher 
(M k = —13) cutoffs in Table [TT] 

In addition to comparisons with observations, we can 
also compute the density of GCs using cosmological sim¬ 
ulations. The publication of the GC Systems catalog in 
m noted a correlation between the dynamical mass of a 
galaxy and the size of its GC population. This relation¬ 
ship was expanded upon in [22], which measured a very 
strong correlation between the mass of the GC popula¬ 
tion and the galaxy halo mass. This relation takes the 
following form: 


logic m gcs =a + /3(log 10 M h - (log 10 M h )) , (5) 

where a is 7.706(7.405)(7.157), /? is 1.03(0.96)(1.21), and 
(log Mh) is 12.3(12.2)(12.2) for all(low-metallicity)(high- 
metallicity) GCs. Unlike the Ncc — relationship, 
the relationship between halo mass and Mqcs does not 
strongly depend on the galaxy morphology. 

In order to convert this to a spatial density of GCs, we 


Mfc < -17 M k < -15 M k < -13 


All 

0.62j;i 

0.77j:i 

0.998:1s 

Elliptical 

n i 7O.32 
U- 1 1 0.09 

n 1 7O.32 
U- 1 '0.09 

o.i78:8l 

Lenticular 

fl 1 o0.25 
u * i ' : >0.07 

A -1 o0.25 
u - lc >0.07 

n 1 9O.25 
u - 1 ' ) 0.07 

Spirals (Sa-Sd) 

o.i38:g! 

0.138:81 

0.138:8! 

Irregular 

o.iiBil 

0.27 8:2? 

0.498:1! 


TABLE II. The number density per Mpc 3 of GCs in the lo¬ 
cal universe, found by combining the Harris GCS catalog m 
with the K-band Schechter Functions from the GAMA survey 
}22j . The errors are found by integrating the lcr uncertain¬ 
ties by the Schechter functions. Note that these errors are 
incomplete, as they ignore the uncertainties in the fits of the 
Schechter functions themselves. 


can multiply this relationship by the dark matter halo 
mass function, as determined by recent cosmological sim¬ 
ulations. We use the functional fit to c u^j h from m , as 
calculated at redshift z = 0 by the HMFcalc website [38] . 
We then compute the integral 


PGC 


M G cs(M h ) dn 

(Mac) dM h 


dM h , 


( 6 ) 


where we use (Mac) = 3 x 1O 5 M 0 , the mean of the GC 
mass function from [26] used in the main text, to convert 
from the mass of a GC system to the number of GCs. 
This yields a spatial density of poc = 3.42 /i 4 Mpc -3 , 
or pgc = 0.72 Mpc -3 , assuming the value of h = 0.677 
used throughout the text. We can also use the similar 
values for low and high-metallicitiy GCs quoted below 
Eqn. [HJ This yields estimates of p l Qc = 0-44 Mpc -3 and 
Pgc 1 = 0-34 Mpc -3 , respectively. 


Models 

This letter considers the BBH inspirals from 48 sep¬ 
arate GC models generated with our orbit-averaged 
Henon-type Monte Carlo code, CMC. The majority of 
these models were first developed in [5]. The full details 
of CMC can be found in previous papers QT[ [321 HO] , but 
the features most relevant to this letter are as follows: 

• three-body binary formation, which we implement 
with a probabilistic analytic prescription [5] > 

• strong single-binary and binary-binary stellar en¬ 
counters, implemented with the small -N integrator 
Fewbody ED, and 

• single and binary stellar evolution with the SSE and 
BSE packages [12] 02]. Note that our implemen¬ 
tation includes several improvements [44], includ¬ 
ing the stellar remnant prescription and BH kick 
physics from [3]. For BBHs which merge within 
the cluster, the GW timescale is calculated by BSE. 
For ejected binaries, the inspiral time is found by 
integrating the orbit-averaged Peters equations [¥5] 
using the masses, separation, and eccentricity of the 
binary at the time of ejection. 

The models begin with 2 x 10 5 , 8 x 10 5 , and 1.6 x 10 6 
number of particles, and are evolved to an age of 12 Gyr 
each. We do not include GC models which dissipate be¬ 
fore 12 Gyr, as we have no way to compare these models 
to observations. However, since these models all begin 
with low numbers of particles and produce low numbers 
of BBHs, the effect on the rate estimate will be minimal. 
For a complete list of GC initial conditions considered, 
see Table [Till 









We also explore the space of initial cluster sizes and 
concentrations, with initial virial radii of 0.5, 1, 1.5, 2, 
and 4 pc and initial King concentrations (Wo) of 2, 5, 
7, and 11. We qualitatively find that the initial King 
concentration does not have a strong influence on the 
observable GC properties at 12 Gyr. Each of our mod¬ 
els starts with 10% of the objects in primordial binaries, 
and stellar masses chosen from a universal initial mass 
function (IMF) [IBl . 

We also explore metallicities of Z = 0.005, Z = 0.001, 
and Z = 0.0005, which are placed at galactocentric dis¬ 
tance of 2, 8 and 20 kpc respectively. This is due to the 
observed correlation between GC metallicity and galacto¬ 
centric distance [47] . Although we explore three distinct 
metallicities, we separate our models into “low metallici- 
tiy” (those for which [Fe/H] < —0.8, i.e. Z = 0.0005 and 
Z = 0.001), and “high metallicity” (Z = 0.005) GCs. 
This is chosen to simplify the comparison between our 
models and the observations of GCs, which show a strong 
bi-modality in metallicity [U]. We also assume that the 
fraction of low-metallicity GCs is 0.76, since that is the 
fraction of MWGCs for which [Fe/H] < -0.8. 
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Initial Conditions 

Properties (12 Gyr) 

N (xlO 5 ) 

R v 

(pc) Metallicity (z) 

Mass (1O 5 M 0 ) 

Rc/Rh 

BH Retained N \ nsp 

2 

0.5 

0.001 

0.55 

0.17 

0.05 

16 

2 

1.0 

0.001 

0.20 

0.34 

0.03 

11 

2 

1.0 

0.0005 

0.68 

0.41 

0.21 

8 

2 

1.5 

0.001 

0.54 

0.42 

0.11 

8 

2 

1.5 

0.0005 

0.62 

0.45 

0.18 

7 

2 

2.0 

0.0005 

0.63 

0.41 

0.20 

5 

2 

2.0 

0.001 

0.55 

0.77 

0.16 

1 

2 

2.0 

0.0005 

0.64 

0.53 

0.20 

8 

2 

2.0 

0.001 

0.56 

1.12 

0.18 

6 

2 

2.0 

0.0005 

0.65 

0.45 

0.23 

5 

2 

2.0 

0.001 

0.56 

0.47 

0.20 

5 

2 

2.0 

0.001 

0.50 

0.99 

0.22 

3 

2 

4.0 

0.001 

0.68 

0.58 

0.37 

2 

8 

0.5 

0.0005 

2.58 

0.27 

0.10 

145 

8 

1.0 

0.001 

2.17 

0.46 

0.18 

91 

8 

1.0 

0.0005 

2.78 

0.48 

0.24 

87 

8 

1.0 

0.005 

2.00 

0.45 

0.14 

88 

8 

1.5 

0.0005 

2.73 

0.47 

0.32 

81 

8 

1.5 

0.001 

2.61 

0.54 

0.30 

67 

8 

1.5 

0.005 

2.02 

0.64 

0.23 

59 

8 

2.0 

0.0005 

2.76 

0.52 

0.43 

48 

8 

2.0 

0.001 

2.62 

0.58 

0.40 

53 

8 

2.0 

0.005 

2.04 

0.63 

0.30 

45 

8 

2.0 

0.0005 

2.79 

0.59 

0.40 

49 

8 

2.0 

0.005 

2.00 

0.52 

0.32 

45 

8 

2.0 

0.001 

2.66 

0.41 

0.36 

59 

8 

2.0 

0.0005 

2.81 

0.77 

0.44 

66 

8 

2.0 

0.001 

2.70 

0.5 

0.40 

57 

8 

2.0 

0.005 

1.92 

0.77 

0.32 

52 

8 

2.0 

0.001 

2.60 

0.67 

0.45 

38 

8 

4.0 

0.001 

2.91 

0.61 

0.58 

27 

16 

1.0 

0.0005 

5.44 

0.5 

0.28 

324 

16 

1.0 

0.001 

4.69 

0.65 

0.27 

269 

16 

1.0 

0.005 

4.49 

0.52 

0.24 

270 

16 

1.5 

0.0005 

5.59 

0.72 

0.42 

204 

16 

1.5 

0.001 

5.39 

0.71 

0.41 

175 

16 

1.5 

0.005 

4.75 

0.59 

0.33 

203 

16 

2.0 

0.0005 

5.68 

0.57 

0.56 

159 

16 

2.0 

0.001 

5.50 

0.56 

0.49 

154 

16 

2.0 

0.005 

4.84 

0.75 

0.43 

162 

16 

2.0 

0.0005 

5.76 

0.6 

0.52 

172 

16 

2.0 

0.001 

5.56 

0.66 

0.52 

139 

16 

2.0 

0.005 

4.97 

0.57 

0.41 

175 

16 

2.0 

0.0005 

5.76 

0.56 

0.53 

167 

16 

2.0 

0.001 

5.58 

0.67 

0.52 

140 

16 

2.0 

0.005 

4.96 

0.63 

0.41 

153 

16 

2.0 

0.001 

5.47 

0.61 

0.52 

121 

16 

4.0 

0.001 

5.94 

0.79 

0.67 

99 


TABLE III. List of the 48 GC models used in this study. The initial conditions are varied across the number of initial 
particles (N), the initial virial radius and the initial metallicities. These models also explore a number of different initial King 
concentrations (u>o), but those are excluded from this table, as they are not observed to have a significant correlation with 
the observed properties at 12 Gyr. We also include the observational properties after 12 Gyr of evolution, including the final 
GC mass, the ratio of the core radius to the half-light radius, the fraction of total BHs remaining in the cluster, and the total 
number of BBHs formed by each cluster that inspiral within 12 Gyr. 
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Advanced LIGO [Phys. Rev. Lett. 115, 051101] 


Carl L. Rodriguez, 1 Meagan Morscher, 1 Bharath Pattabiraman, 1,2 
Sourav Chatterjee, 1 Carl-Johan Haster, 1,3 and Frederic A. Rasio 1 4 5 

1 Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) and Dept, of Physics and Astronomy, 
Northwestern University, 21f5 Sheridan Rd, Evanston, IL 60208, USA 
2 Dept, of Electrical Engineering and Computer Science, Northwestern University, Evanston, IL, USA 
3 School of Physics and Astronomy, University of Birmingham, Birmingham, B15 2TT, United Kingdom 

(Dated: December 21, 2015) 


We have discovered an error in our integral over the globular cluster mass function that was reported in the 
letter. The computational error artificially over-counted the contribution from high-mass globular clusters, artificially 
increasing the mean number of mergers per globular cluster over 12 Gyr. This changed the quoted values in Table 1, 
and also increased the quoted detection rate in the paper. The corrected Table 1 is reproduced below 



Mass Cutoff 

Metallicity 

4 x 1O 6 AL 0 

2 x 1O 7 M 0 

2 x 1O 8 M 0 

Low 

159 

223 

250 

High 

238 

424 

549 


TABLE I. The corrected mean number of inspirals per GC over 12 Gyr. 

With the updated values from Table 1, the quoted detection rate for Advanced LIGO decreases. The realistic rate 
of detections is now 30 mergers per year. The pessimistic rate is 10 mergers per year. The optimistic rate is 100 
mergers per year. 

Fortunately, this correction significantly decreases the uncertainties quoted in our original letter. Originally, the 
dominant uncertainty was the upper-mass cutoff of the globular cluster mass function; however, by correctly accounting 
for the contribution of high-mass clusters, the upper-mass cutoff no longer heavily influences our rate estimate, and its 
associated uncertainty is significantly reduced. The dominant uncertainty in the updated estimate is now the spatial 
density of globular clusters per Mpc 3 in the local universe (see the first Supplemental Material). As noted in the 
letter, the spatial density used in previous estimates of the binary black hole merger rate is ~ 2.4Mpc -3 [1], similar to 
the number we employ in our optimistic rate estimate. As such, the rate we find from our models, when normalized 
to the same value of pec , is still ~ 5 times that of previous estimates [2]. 

This correction does not change the major conclusions of the letter, that dynamically-formed binaries contribute 
significantly to the overall merger rate, and that those binaries are characteristically more massive than those formed 
by pure stellar evolution. We reiterate that the mass difference between binaries formed dynamically and those formed 
in the field is heavily dependent on the physics of binary stellar evolution. In particular, we employed the original 
stellar wind prescription of [3], whereas more recent studies (e.g. [4]) have used the enhanced stellar wind prescription 
of [5]. Efforts to understand the contributions of varying binary stellar evolution physics are currently underway. 

We are very grateful to Cole Miller for pointing out this error. 
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